Non-Equilibrium Statistical Physics of Currents in Queuing Networks 



Vladimir Y. Chernyak "'''j^] Michael Chertkov "'''jT] 
David A. Goldberg "'"^j^j and Konstantin Turitsyn "''j^ 

Center for Nonlinear Studies and Theoretical Division, LANL, Los Alamos, NM 87545 
''Department of Chemistry, Wayne State University, 5101 Cass Ave, Detroit, MI 48202 
'^New Mexico Consortium, Los Alamos, NM 87544 
Operations Research Center, MIT, Cambridge, MA 02139 
Landau Institute for Theoretical Physics, Moscow, Russia, 119334 
(Dated: June 22, 2010) 

We consider a stable open queuing network as a steady non-equilibrium system of interacting 
particles. The network is completely specified by its underlying graphical structure, type of inter- 
action at each node, and the Markovian transition rates between nodes. For such systems, we ask 
the question "What is the most likely way for large currents to accumulate over time in a network 
?" , where time is large compared to the system correlation time scale. We identify two interesting 
regimes. In the first regime, in which the accumulation of currents over time exceeds the expected 
value by a small to moderate amount (moderate large deviation), we find that the large-deviation 
distribution of currents is universal (independent of the interaction details), and there is no long- 
time and averaged over time accumulation of particles (condensation) at any nodes. In the second 
regime, in which the accumulation of currents over time exceeds the expected value by a large 
amount (severe large deviation), we find that the large-deviation current distribution is sensitive 
to interaction details, and there is a long-time accumulation of particles (condensation) at some 
nodes. The transition between the two regimes can be described as a dynamical second order phase 
transition. We illustrate these ideas using the simple, yet non-trivial, example of a single node with 
feedback. 

Keywords: Statistics of Non-Equilibrium Currents, Open Queueing Networks, Condensation phenomenon, 
Birth-Death Processes 



I. INTRODUCTION 



A. Non-equilibrium statistical physics and queueing networks 



The concept of statistical equilibrium is extremely powerful. Once detailed balance (which is 
synonymic to the equilibrium) is established, one can shortcut a discussion of dynamics and just 
consider the Gibbs distribution that governs simultaneous correlations in the steady state. On the 
other hand, if detailed balance is broken no free lunch is guaranteed, and one generally does need 
to dive into dynamics, even to describe just the steady state. This is an infamous and principal 
difficulty at the very core of non-equilibrium statistical physics. It is thus of interest to identify a 
class of non-equilibrium steady systems where the steady state, distinctly different from the Gibbs 
distribution, can be derived in a straightforward way. 

The problem discussed in this manuscript belongs to this class - it is a non-equilibrium statistical 
physics problem with the steady state known explicitly. More precisely, we study the model generally 
known as an open queueing network. A general (Markovian) open queueing network, the so-called 
Jackson network [T], can be described as a random walk of particles (jobs, vehicles, people, computer 
packets, etc) on a directed graph. The directed links label the transitions between stations/nodes, 
each of the Poissonian type and thus characterized by a single number (the rate). Each node is 
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FIG. 1: Example of an open finite queueing network represented by a directed graph. The sample graph consists of four 
vortexes/stations, labeled 1,2,3,4, with label is reserved for an external (out) node. Transitions between the stations are 
shown as directed edges. Loops (self- loops), as 1 — > 1, are allowed. Each graph edge is equipped with a transition rate. 
Throughout the manuscript all transitions in the network are assumed memoryless, i.e. 0-degree Markovian - Poisson. We also 
focus on the case of the infinite waiting room, i.e., no particles/jobs are lost, and thus the number of particles accumulated at 
a node may reach +00. Another important characteristic of a station/note is the number of tellers/servers. For the sake of 
simplicity, we assume the jobs/particles to be identical (single-class), and all tellers at any given station processing with the 
same time-independent Poisson rate. 



characterized by the number of equivalent servers. Collectively this defines a many-particle prob- 
lem, whose non-equilibrium nature (no detailed balance) follows immediately from the definitions. 
Adopting the terminology of the Queueing Theory community (a subset of the Operations Research 
community), the type of service at a node is provided by the M/M/m/oo queue, which is translated 
as Markovian input, Markovian output with m servers, and infinite waiting room. (Throughout the 
manuscript we use the shorter notation M/M/m.) An example of such a network is shown in Fig. [l| 
and more information will be provided below. Note that a similar, and to a degree more general, 
model is known in statistical physics as the zero range model [2J. 

In spite of the generally pessimistic non-equilibrium assessment, the stable (i.e. achieving a statis- 
tical steady state) open Jackson network allows an explicit and simple solution for the steady state 
[3]. Similar statements apply to the zero range model |2j. The solution for the steady state is the 
so-called product form [2]-S] , where the joint distribution function of the occupation numbers at all 
the nodes/stations is factorized into a product of marginal probabilities at the separate nodes. Some 
new and physics-based exposition of this factorized solution is due to and will also be a part of 
our construction below. We note that a different factorization has also led to the solution of the 
related Asymmetric Exclusion Processes (AEP) Models, see P-[9] and the references therein ^ . 

On the other hand, in recent years the quest for universality in non-equilibrium statistical problems 
has turned to the analysis of currents generated over time, long compared to the system correlation 
time scale. This route became fruitful and helped to establish some fundamental relations about 
the symmetry of the currents distribution, known as fluctuation theorems [T014T2] . In spite of this 
partial success, the full description of the currents distribution for a general many-particle and open 



The general AEP models are significantly different from the Jackson-network models discussed in this manuscript - they can be viewed 
as a special case of an M/M/1/1 queuing model (waiting room with one slot and no particles lost) in contrast with the Jackon network 
of M/M/m/oo (infinite waiting room) queues. In the case of a simple one-dimensional chain the AEP model is reducible (with a proper 
redefinition of the phase space) to the Jackson network model (and other way around) [6]{9], and then both models show product-state 
form solution [2}jl]. However, the product-state solvability, which holds for the Jackson network generally, does not extend to the AEP 
model on general graphs. 
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non-equilibrium problem was (and is) deemed too difficult ^. 

B. Large deviations, queueing networks, and condensation phenomenon 

One of the central questions in the Queueing Theory community is the probability of rare events in 
queueing networks (large deviations) [15]. A central tool in understanding such events is the notion 
of the so-called 'fluid-limit' of a queueing network [TB], which was originally developed to understand 
certain questions related to the stability (existence of steady state) for certain complicated networks 
[T71 ITS] . We now describe this fluid-limit in greater detail. For a fixed queueing network, one scales 
both time and the initial number of particles in system by some large integer n, and then 
normalizes the set of queue lengths by n. Note the important difference with the standard notion 
of hydrodynamic limit - here, the number of nodes in the network is fixed - only time and the 
initial number of particles are scaled. It is proven in [151 IB] that for a broad class of networks, 
this scaling has a non-trivial limit (the fluid limit). It is then standard in the Queueing Theory 
literature to perform the large deviations analysis in the setting of this fluid limit, as opposed to 
the original (unsealed) network. There has been great progress in understanding all aspects of 
these large deviations [T9H23] . However, driven by applications to describing buffer overloads in 
communications systems, most of these results have been geared towards understanding how large 
queue lengths accumulate over time (see for example [22] )• Much less is known about how large 
currents build up in a queueing network over time, especially for the unsealed (not fluid limit) 
network. 

The statistical physics community has developed several tools ( e.g. the fluctuation theorems 
[TUHT^ ) that are well-suited to studying the large-deviations properties of currents in a network. In 
light of the aforementioned gap in our understanding of the large deviations of currents in queueing 
networks and related network models, several researchers in the statistical physics community have 
recently begun to apply the fluctuation theorems to the study of currents in a variety of network- 
related models [231 - I27] . Most of these results have been for models such as the zero-range process, 
which are tangential (although closely related) to the networks studied by the Queueing Theory 
community. We bridge this gap by directly applying such an analysis to the canonical model of 
queueing theory - namely the Jackson network. We note that although the potential application 
of these tools to Jackson networks is mentioned in [21], the paper closest in spirit to our own 
is [25]. Indeed, for the 1-D zero-range process the authors identify several regimes in which the 
model may operate, characterized by the large-deviations properties of currents, and whether or 
not there is an accumulation of particles at sites over time (so-called condensation phenomenon). 
Their analysis proceeds by formalizing the system dynamics in a quantum-mechanical / operator- 
theoretic framework, and then studying the relevant spectral properties and Cramer (large deviation) 
functions. Our own analysis will very-much parallel that of [25], but for the Jackson network model 
on general graphs. Other related work can be found in [2H], in which the zero-range model on regular 
lattices is studied in the hydrodynamic limit, as well as [22] where the closely related stochastic lattice 
gas model is studied. 

We also note that the Queueing Theory community has already had some success in applying 



^ We would like to emphasize here the principal difficulty arising from the fact that the system is both (a) many-particle and (b) open. 
The analysis of such systems is often limited to the setting in which there are few degrees of freedom. An example is a recent paper of 
three of us with Puliafito 13: that discusses an explicit expression for the distribution of currents associated with a polymer stretched by 
external shear flow. In another paper (two of us with Malinin and Teodorescu 114J) considered a setting with multiple degrees of freedom 
and analyzed the statistics of so-called topological currents. The statistics of currents have also been studied for the aforementioned 
AEP models over one-dimensional chains, which typically have multiple degrees of freedom [S]- 
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some of these tools [301432] . Indeed, these analyses (and our own) rely on the interpretation of 
queueing systems as interacting particle systems, which (historically) helped lead to many of the 
great breakthroughs in the understanding of queueing networks (e.g. product-form solution [33]). 
However, these analyses have been for either an infinite 1-D chain of queues [301 E2] or only in the 
hydrodynamic limit [HI]. In either case, their findings parallel our own, in the existence of certain 
phase transitions in which condensation may (or may not) occur. 

C. Currents in a queueing networks with feedback 

We now discuss the nature of the currents in a Jackson network in steady state, a topic that has 
generated much research in the queueing literature [3H - I38] . Although currents over the entering links 
are Poissonian by construction, those which leave the system are not obviously Poissonian. However, 
if the system is stable, these exiting flows are in fact Poissonian [33l [35] . More generally, the flow 
is Poissonian along all arcs that may not be revisited by a particle, and these flows enjoy several 
nice properties (such as asymptotic independence) [38]. It was recognized early on in the queueing 
literature that the statistics of the internal currents with feedback are rather difficult to analyze [33] . 
The source of these difficulties is the complicated feedback mechanism that arises when particles may 
revisit an arc. Thus, even for the simple single-node feedback system (which will serve as an enabling 
example in this manuscript), the statistics of the feedback current are rather complicated [39tiH] . 
Although several properties of the currents in a multi-server feedback queue and its generalizations 
have been studied in the literature [ISHU], most of these results involve showing that in certain 
limiting regimes the flows are close to Poissonian under some metric [131 H8H50] - the precise nature 
of these flows remains poorly understood. Furthermore, it seems that an operator-theory/generating- 
function approach has not yet been applied to the study of these feedback currents. We will take 
this approach to derive new understanding of these currents. 

D. Our main results 

The main results of our manuscript can be described as follows: 

• We present an explicit, detailed operator-theoretic description of a Jackson network in the Doi- 
Peliti formalism. This expands on the description given in [5], with an eye towards introducing 
a larger set of the statistical physics community to standard queueing models. 

• We identify an "uncongested" regime in Jackson networks, w.r.t. the large-deviation behavior 
of current. This regime is universal, i.e. interaction-independent, and non-Poissonian. Further- 
more, in this regime there is no infinite accumulation of particles at any node (condensation). 
This universality can be qualitatively explained as follows. In this regime, the given deviation 
is driven by sample paths in which the number of particles does not diverge. This will generally 
occur when the given deviations are somewhat mild, and a deviation can be attained without 
a massive buildup of particles. Then, the impact of one particle 'blocking' another does not 
contribute asymptotically to the deviation, and thus the system is equivalent asymptotically to 
one in which all nodes have an infinite number of servers (particles do not interact). We confirm 
the existence of this regime by demonstrating that the single-node network with feedback falls 
into this regime for certain parameters, which we compute explicitly. 

• We also identify a second, "congested" regime, which is interaction-dependent. In this regime, 
a given deviation is driven by system primitives that are dependent on the number of servers 
and service times. This will generally occur when the given deviations are more severe, and the 
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only way to attain the given deviation is to have all, or at least some, servers busy for essentially 
the entire time horizon. In this regime the time-averaged queue length diverges for at least one 
node, marking the dynamical phase transition between the two regimes as second-order. 

• We observe that statistics of time-averaged (over the large observational interval) queue and 
queue measured at the last moment of time, both conditioned to an atypical (large or small) 
values of currents, are not identical. The difference is particularly striking in the "congested" 
regime of the largest currents, where all moments of the former object (queue at the last 
moment of time) saturate to finite values while all moments of the later object (time-averaged 
queue) diverge with time. Quite generally this phenomenon can be classified as a breakdown 
of ergodicity in cases which are atypical with respect to currents. 



E. Outline 



The manuscript is organized as follows. In Section [TT| we present a technical introduction to the 
dynamics of the Jackson network in terms of the physics-native Doi-Peliti ("quantum" or "second- 
quantized") technique [5TH53] . There we formally characterize the Master Equation (ME) that 
governs the system evolution and steady state, as well as the joint distribution of densities (occu- 
pation numbers). We customize this description to the M/M/oo, M/M/1, and general M/M/m 
models in Sections II A , II B and II C| respectively. In Section [ITU which represents the core of the 
manuscript, we adopt the Doi-Peliti technique to analyze the ME for the joint distribution function 
of densities (that reside at the nodes) and currents (that reside at the links). We also show how the 
coherent-state technique provides a complete description of the ground state (eigen- value and eigen- 
function) in terms of the relevant evolution operator. Section III is partitioned into three Subsections 
and an Appendix. In Section III B , we discuss the universal (and statistically typical) "uncongested" 
regime. We also describe the boundary of the "uncongested" region in the space of currents, and 
comment on the associated dynamical phase transition. Our analysis procedes by invoking an auxil- 
iary construction for the left eigen-function of the evolution operator, which is discussed in Appendix 
[a} The transitions between the "uncongested" and (partially) "congested" regimes are discussed in 
Section III C The general theory is illustrated in Section HID for a single-node system with feed- 
back. In Section IV, to validate the theory, we describe a full spectral solution for this single- feedback 
problem. In Section W\ we draw several conclusions and discuss future directions for research. 



II. THE DOI-PELITI-MASSEY OPERATOR TECHNIQUE FOR A GENERIC BIRTH-DEATH 

PROCESS 



This Section introduces notation and describes the main operational rule of the Jackson network 
in terms of the statistical physics native Doi-Peliti technique We note that a very similar 

(but less explicit) formulation was derived in [5]. 

As we will see, in this context the Doi-Peliti technique is closely related to the operator-theoretic 
framework formulated by Massey [5^ [55] in the Queueing Theory community. We start by introduc- 
ing the quantum-mechanics based bra(c)ket notation. We then discuss the product-form solutions 
for the stationary problems associated with the M/M/oo, M/M/1, and generic M/M/m networks 
in Sections II ApiB and II C[ respectively. 

The network ( e.g. the one shown in Fig. ([T]) ) is represented by the directed graph, [Qq^ ^i), where 
Qq,Qi marks the set of vertices and directed edges of the graph (respectively). If at some instance t, 
node i has a queue of size n, one says that the node is in the state represented by the ket- vector 
where n = 0, 1, ■ ■ • . Then, any "pure" state of the network will be denoted by the ket-vector |n). 
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where the components rii of a vector n = {ni\i G Qo) are labeled by the network nodes (vertices). If 
a state \n) is realized with the probability P{ti), we say that the entire network is in the following 
"mixed" state 

|.) = 5^P(n)|n), $^PH = 1, (1) 

n n 

where the last condition reflects the fact that the total probability equals unity. Here and below we 
formally assume that P{n) = whenever any component of the vector n is negative. 

It is convenient to introduce a Hilbert space of ^o-dimensional analytic functions of the vector 
variable z = {zi\i G ^o) 

P(z) = 5^P(n)n^r, (2) 

which is also known as the generating function in the theory of birth-death processes |56l EZ] . 
The "quantum" (pure) states are transformed by the following creation and annihilation operators: 

a+l ■ ■ ■ ■ ■ ■ ) = I ■ ■ ■ + 1, ■ ■ ■ ), aj\--- , n^, ■ ■ ■ ) = | ■ ■ ■ ,nj - 1, ■ ■ ■ ). (3) 

The normalization condition in Eq. ([T]), i.e. conservation of probability, reads 



(0|exp I Xl^i I 1^) = 1' 



(4) 



where the vacuum state |0) = |0, ■ ■ ■ ,0) corresponds to the empty queue over the entire network. 
In these notations ME becomes 

dt\s) = H\s), (5) 



where H is the Hamiltonian operator of the Q-network. In an integrated form, Eq. (|5j) is equivalent 
to 



\s{t)) = U{t)\s{0)), U{t) = Texp (^j^ dt'H^ , (6) 



where Texp is defined as a time-ordered exponential, i.e. the product of time-discretized operators. 
Furthermore, it becomes normal exponential if the parameters of H (i.e. transition rates) do not 
carry explicit time dependence. 

Note that the Hamiltonian H is always real, as it represents probabilities which are positive and 
bounded. Thus the normalization conditions Q and Eq. ([s]), which should be enforced by the theory 
for any feasible n, result in 



(7) 



where we have used standard bra-vector notations. Recall that an operator acting on the bra- vector 
from the right generates a bra-vector, and all the features of left operations can be extracted directly 
from the normal definition of the bra-(c)-ket scalar product, ^n,m : {n\m) = 5{n,m). Stating it 

differently, (0| exp \ % ) is the left eigen- vector of the Hamiltonian with zero eigen- value. 



The expectation value of a (dummy) operator • over a state \s) is 



») = (0| exp I XI I 



and according to Eq. ([T]), the corresponding Heisenberg (evolution) equation becomes 

d,{i) = {[i,H]). (9) 

Here we have assumed that • does not have an explicit time dependence, and [A, B] is the standard 
notation for a commutator. 

A. M/M/oo network 



The generic form of the ME for a network of M/M/oo queues is 

— P(n; i) = X {{ui + 1)P(- ■ ■ ,ni + !,■■■ , - 1, ■ ■ ■ ;t) - niP{- ■■ ,ni,--- ,nj,--- ;t)) 
+ J2^oi (P(-- - ,n,-l,--- ;t)-P(... ;t)) 



+ J] A,o ((n^ + 1)P(- • ■ ,ni + l,--- ■t)-n,P{--- ,n„----t)). (10) 



Here stands for the directed edge of the network corresponding to a job transfer from site i to 
site j, with Poisson rate Xij; and Xoj,Xjo are the Poisson rates for job injection and removal to/from 
the network at site j (respectively). Applying summation over properly weighted n-states to both 
sides of Eq. (10), and using the relations Eq. (5][6), one arrives at the following Hamiltonian: 



Hoo= ^ Aij(a+ - a+)ai + X (Aoi(a+ - 1) + Aio(l - a+)ai) , (11) 
(j,i)eei i6Go 

which was first derived for the problem in |52j . 

We further introduce a path-integral representation. The analytic structure of the theory is as 
follows: 

t) = j ^^W{z, C)V{C; 0) exp(-CC'), (12) 

r'n'(t)=z / /■* \ 

W{zX)= VrjVrj'expizriit)- dt' [r^' {t')r){t') - HM {t')Mt)))] . (13) 

Jvio)=C \ Jo J 

where Tiooiil'^ V) corresponds to Hoc expressed as a polynomial over the creation/annihilation oper- 
ators, so that the creation operators are all positioned on the left from the annihilation operators 
(normal ordering), and a^,aj are replaced by rj'^ and rjj (respectively). Thus, for a general birth-death 



model (11) over a network ^, one derives 



V-oo{r}\rj)= J2 A,,(r7;-r/:)r7, + 5^(Ao.(77:-l) + A,o(l-r/:)r/,). (14) 



As usual, the path-integrals in Eq. (13) should be understood as the continuous limit of the following 



discretized multiple integral (see [52j for explanations and accurate validation of the proper discrete- 
time regularization) : 



„ Af-l 

W{zX)=\im / n 

^ 7 1 



dr]idr][ 



X 



J- (27rz)l^?ol 

exp [ z'nN-i + AH(z, r]N^i) + ^ {-r]'i{r]i - m-i) + ^'^i'n'i, m-i)) 



(15) 



where A = t/N. 

Finally, one finds that the creation-annihilation (11) and the path-integral ( 13p5 ) formulations 
of the birth-death process also allow for the following simple "differential" interpretation for the 
analytic function defined in ([2]): 

dtViz;t) = nM'Piz;t), (16) 

iioo{,z) = nooiz, 9^) = ^ Xij{zj - z,i)dz^ + (Xoiizi - 1) + Aio(l - Zi)dz^) . (17) 

(i,i)egi i&Go 

Thus the mapping from the creation-annihilation operators to the poly-differential operators (holo- 
morphic representation) is 

(a+,a) ^ {z,d,). (18) 



Looking for a steady state (time-independent) solution of Eq. (17) in the exponential form 



Poo (2) = exp 



hi{zi - 1) 



and substituting the ansatz into Eq. (16), one arrives at the following set of conditions on h: 



(Aoi - Kohi) = 0, 

i&Qo 

Vi G ^0 : - hi Xij + ^ Xjihj + Xoi - Xiohi = 0. 



(19) 



(20) 



(21) 



Although it seems that there is one more condition than the number of variables, the conditions 



are dependent (sum Eq. (21) over all vertices of the graph). Therefore, solving the system of 
inhomogeneous linear equations Eq. (21), which we restate for convenience as 



A/l ^in 



-XoilieQo), A = (Aij|i, J G ^0), Aj 



Xji, 



consists in evaluating 

h = A-'Xin. (23) 
Here the existence of the steady state solution requires that: (a) A is not singular, and (b) all 



components of the h vector that solves Eq. (23) are positive. 
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Note that the form of Eq. (19) is fully factorized. Thus once the valid solution of Eq. (21) is 



found, the full probability of observing the system in any given state is decomposed into a product 
of probabilities, each evaluated at the relevant graph vertex. Recall that this occurs in spite of the 
fact that to find the re-normalized rates hi one must solve a graph-global linear problem. This strong 
symmetry of the Poisson-In-Poisson-Out process, observed in spite of the fact that the DB is broken, 
is referred to (in the Queuing Theory literature) as "quasi"-DB (TS^, ^33]. 

The special feature (memoryless property) of the exponential distribution is also very transparent in 
the creation-annihilation language. Indeed, one observes that the exponential (in quantum mechanics 
also referred to as "coherent" ) state exp(/ia+)|0) is the eigen-f unction of the annihilation operator 
a, with the eigen-value h 



d|cohoo(/i)) = h\cohoo{h)) , \cohoo{h)) = exp{ha'^)\0) . 



(24) 



Therefore, 



^oolcohoo(^)) = I Yl ^'^^^t - + (Aoi(a+ - 1) + Aio(l - a+)/i,) | |cohoo(^)), (25) 



and the stationarity condition Hoo\co\ioo{h)) = translates exactly into Eqs. (21), where the i-th 



equation correspond to the condition that the c-factor in front of the corresponding dj' is zero. 



B. M/M/1 network 

Consider a network of M/M/1 processes. In this case the ME adopts the following form: 
d 



dt 



P{n; t) = Y A,; J P(- ■ ■ , rii -Fl, ■ ■ ■ , rij - 1, ■ ■ ■ ; t) - P(- ■ ■ , nj, ■ ■ ■ , rij, ■ ■ ■ ; t) 



(i,i)6gi 



t)) 



+ ^Aio (P(--- ,ni + l,--- ■,t)-P{--- ,ni,--- ■,t)). 



(26) 



Here d{x) is the characteristic function of the logical condition x, i.e. it is unity when the condition 
is satisfied and zero otherwise. The corresponding Hamiltonian operator in Eq. (5][6) is of the form 



Hi= Y^ '^^jK^'j 

(ij)egi 



at — at 



(27) 



Here hi is a "skewed" annihilation operator (see e.g. [5] for a similar operational rule), such that 
bi\ni) = 9{ni > 0)|nj — 1), and in both Eq. (26) and Eq. (27) we keep the same notation as in 
Eqs. ( lOpT ) (respectively). 

Note that b is expressed in terms of a and d^ in an extremely nonlinear way. However, the 
representation allows a simple "analytic" interpretation [5]: 



^^^Pn|^) — ^ ^-^-^ — P^S})_^ where p{z) = 



(2^ 



10 



For the introduced generating function representation, the analog of Eqs. (T6|T7) becomes 



dtV{z;t)= Xij{zj-Zi 



V{z-t)-V{z^i]t) 



+ $^A,o(l 



- Zi 



V{z;t)-V{z^f,t) 



(29) 



Here z^j = ((1 — 6ij)zj\j G ^o)- In words, this is the vector z with the component Zi replaced by 



zero. 



Let us come back to the skewed-creation-annihilation representation, and note that the coherent 
states associated with this "skewed" annihilation operator b were discussed in |5j . The approach can 
also be traced back in the Queuing Theory literature to the classic papers of Massey [Ml ES] on the 
operator approach to Jackson networks. The coherent states for b are constructed as follows: 



6|cohi(/i)) = /i|cohi(/i)), |cohi(/i)) 



1 



1 — ha' 



-|0). 



(30) 



Then the analog of Eq. (25) becomes 



ii'i|cohi(^)) 




ij[a. - a, 



-)h, + J2 i^o^{at - 1) + A,o(l - di)h,) I |cohi(/i)). (31) 



Furthermore, the condition of stationarity, Hi\coh.i{h)) = 0, translates exactly into Eqs. (21), where 



the i-th equation corresponds to the condition that the c-factor in front of the respective aj' is zero. 
We conclude that the stationary distribution of the M/M/l-network is 



^i(-) = n I 



1 - hi 



hjZi 



(32) 



where h is the solution of Eqs. (21) 



C. M/M/m network 

The ME in the general case of an inhomogeneous M/M/m- network, with positive integers mi 
(number of tellers) assigned to each vertex i of the graph, can be represented by 

^P{n;t)= J2 0^,(ni + l)^(n, >0)P(--- + 1,--- -I,--- ;t) 

-9mXni)P{- ■ ■ ,ni,--- ,nj,--- ; t)^ 

+ J]Ao. {e{ni>0)P{--- ,ni-l,-- - ;t)-P(-- - ,ni,-- - ;t)) 

+ 5^Aio(^m,(rii + l)P(--- ,n, + l,--- ;t) -0^^(ni)P(--- - ;t)), (33) 

dm{n) = min(n, m). (34) 
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The evolution operator (Hamiltonian) becomes 



H= A.,(a;-a+)r^) + 5^(Ao.(a+-l) + A,o(l-a+)6^ 

fe^'"V) = ^m(n)|n-l). 



{rrii) 



(35) 



(36) 



Thus the problem of finding the stationary solution is reduced (pretty much like before in the m = oo 
and m = 1 cases) to constructing coherent states for the annihilation operator fe^™); 



U'^^\co\im{h)) = h\coKa{h)), |coh„(/i)) = gm{ha-^)\0), 



E 

k=0 



m—1 



X 



m 



k=0 



m 



m—k 



k\ 



(37) 
(38) 



Finally, the full expression for the generating function of the stationary solution over the general 
network becomes 



Viz) 



n 



(39) 



where h is the solution of Eqs. (21). Therefore, by Eq. ([2]), 



(40) 



Obviously, Eqs. ( 39p0 ) are consistent with Eq. (19) and Eq. (32) when m = {mi\i G Qo) is set to 
m = oo and m = 1 (respectively). 

Note (for the sake of accurateness) that when deriving Eq. (40) in the operator formalism we 
took advantage of the important fact that both left (bra-) and right (ket-) zero eigenvalues of the 

Hamiltonian (35), described by (0|exp (^jaj^ and Yli&go 9mi{hia~^)\0) respectively, are explicitly 

known. 

Note that one can recalculate any moment of n from either Eq. (39) or Eq. (40). In particular, for 
the first moment at a station we arrive at 



{n^) 



(0| exp (j2jego «i) ^t^i Uk&Go 9m, {hkd+)\0) q 



9mi {h-i^i 



(0|exp (Eieeo ^i) Y\.k<^Qu 9m,(hka+)\Qi) 



keGo ■ 



dzi gm^hi 



(41) 



z=l 



where gm{x) is taken from Eq. (38) and (as before) h is the solution of Eqs. (21). Note that 
the moments are finite only if hi < rrii, which thus defines the condition for Q-network stability 
(statistical stationarity) [15] . 



III. STATISTICS OF NETWORK CURRENTS 

A (general) queueing network is naturally characterized by the actual current of particles/jobs 
going through and being processed at each station according to a certain service discipline ( e.g. First- 
In-First-Out (FIFO) ) . Here the quasi-current characterizes the activities of tellers, not focusing on 
the dynamics of the individual particles at all. In other words, actual current tracks the dynamics of 
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the jobs/particles, while quasi-current tracks quasi-particles/jobs assuming that all the jobs waiting 
for service at a station are fully equivalent and not prioritized. Actual current and quasi-current 
coincide in the case of an M/M/oo network, when the individual jobs do not interact at all, as well 
as any network in which all particles (customers) are identical. With this disclaimer, we will be 
discussing quasi-currents for the remainder of the paper, and refer to them as currents to simplify 
our exposition. 

As shown below, calculating the statistics of the currents in a Jackson network becomes tractable 
in the two regimes that we will study: the "uncongested" and "congested" regimes. 

In the "uncongested" regime, this arises from the product-form symmetry characterizing the regime 
(an extension of the product-form symmetry discussed earlier), as well as the m- independence 
property (universality). Note that for networks without feedback, the statement of m- independence 
is equivalent to the fact that in steady state, the flow along arcs will be Poisson (with rate independent 
of the number of servers). As mentioned previously, this phenomenon was discovered earlier in the 
Queueing Theory literature ESI ES] • However (and to the best of our knowledge), the extension 
of this statement to the asymptotic (large time) limit for more general networks (unsealed, not in 
the fluid limit or hydrodynamic limit) has not been formally explored in the literature. 

The remainder of this Section is partitioned into four Subsections. We start from a general dis- 
cussion of the current related objects in Section [III A In Section III B we consider the uncongested 
regime. We also develop several generalizations of the Doi-Peliti technique to account for currents, 
and develop some machinery necessary for the statement of our results. At the end of the Section, 
our analysis naturally leads to the identification of the uncongested regime's breakdown, namely 
the identification of a phase transition in the space of currents. We also show that this transition 
is second-order, thus translating into smoothness of the Cramer function of currents (continuity of 
the first and the second derivatives) at the transition. In Section III C , we extend the coherent state 
formalism to the "congested" regime via a simple reduction of the network graph. We note that a 
similar decomposition was applied in [25] , and is similar in spirit to many such reductions appearing 
throughout the Queueing Theory literature [601466] . 

Finally, in Section III D| we illustrate the general theory using our enabling example of a single 
node with feedback. 



A. Preliminary General Remarks 

We will mainly be interested to evaluate the joint distribution function of the currents and queue 
sizes where the latter are averaged over the entire time horizon. In the following we will use P(n, J\t) 
notation for the main object of interest. However it is technically more convenient to start from 
another (and to a degree auxiliary) object defined as a joint distribution function of currents and 
queues where the latter are observed at the final moment of time. We will see below that it is 
important to differentiate these generally distinct objects. 

Let P(n(t), J|t) denote the joint probabihty distribution function of the queue size at the final 
moment of time t, n{t) = {ni(t)\i G ^o); and currents accumulated over the [0;t] interval of time, 
J = (Jijlihj) ^ Gi), where the latter are defined on all edges of the graph and the former are defined 



(as before) on vertices. The ME for this object is the natural generalization of Eq. (33), which we 



now present in operator form (to allow for more compact notations). In particular, the operator 
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form of the ME for P{n{t), J\t) is 



dt\s{n{ty,J)) 




|s(n(t); J)), 



Joi = Aoi(l — a^i)af, Jio = Ko{0'to '~ 1)^: 



(mi) 



(42) 

(43) 
(44) 



where H is defined in Eq. (35). Here we have assumed that the currents are discrete and positive, 



and the respective ket-vector is related to the joint PDF of n(t) and J as follows: 

|s(n(t); J)) =P(n(t);J)|n; J). 



(45) 



Jij ( in Eq. (43) ) is the operator for the amount of current from site i to site j, and afj is the newly 



introduced creation operator ( at edge (i, j) ) acting on the space of discrete positive currents. We 
define operators for incoming and outgoing currents in Eq. (44) similarly. Formal solution of Eq. (48) 
is 



\s{n{t); J)) = exp I t I ^ + ^ 



|s(n(0);J)), (46) 
where the ket-state on the rate correspond to the "initial" steady distribution of queues, described by 



Eq. (40), and zero initial current: |s(n(0); J)) = |s) ® | J = 0), where we follow notations introduced 
in the introduction and \s) = -^(^)l^)- 

It follows that the generating function over the currents € Qi, accounting for the incoming 

{{0,i) E Gi) and outgoing ((i,0) G Qi) arcs, is 



Mnm = J2 n 4^>Ht);J)). 

J (i,i)6Gi 



(47) 



According to our standard birth-death [creation/annihilation] rules, the object described by Eq. (47) 
satisfies 

dt\sq{nit))) = HqMn{t))), \sq{n{t))) =exp{tHq)Y,Pin)\n) =exp{tHq)\s), (48) 



H. 



am 



(j,i)egi {o,j)eei 



1)+ 5^ Aio(g 

{j,o)egi 



iO 



(49) 



Note that even though our evaluation in Eq. (48) applies to the general object, namely to the joint 



distribution function of currents over the entire network, one may compute the relevant marginals 
(say the distribution function for the current over a single edge) in a straightforward manner. We 
also note that the operator Hq can be obtained by modifying/twisting the evolution operator H, 
given by Eq. (p6j), as follows. One weights the off-diagonal terms of H in the space of populations |n), 

namely \ija^of^'\ Xoiaf, and Ajo^,-™''*, with the factors qij, qoi, and qio respectively. This corresponds 
to viewing the underlying stochastic process as a Markov chain on an infinite graph, whose nodes 
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are labeled by pure states |n), whereas links represent the set of processes allowed by the evolution 
operator H. Within such a picture, the set of q parameters plays the role of discrete gauge fields 
(vector potentials), and Hq is interpreted as the evolution operator, "twisted" by the gauge field g, 
as described in [67]. Since here we are dealing with oriented graphs, no constraints are imposed on 

q 

It is also useful to consider the distribution function of queue at the finite moment of time, 
conditioned to specific value of the current generating parameter q. This object, and the respective 
first moment of queue size, become 

Pq{n{t)) oc (n(t)l exp (tH^ (50) 

(0| exp ( Ejeso ) ^^t^i exp (tHA \s) 

Mt)), = \ \ . ^ ^ • (51) 

(0| exp [Y^j^Qo %j [tHq j \s) 

Returning back to our main object of interest (the joint distribution function of the current and 
of the queue averaged over the time horizon) and following the same formalism/notations, we derive 
the analogs of Eqs. ( SOfsI ) 

Pq{n)(xt-^ I rft'(0|exp ( j exp (^(t-tO^g) |n(tO)(n(tO|exp (^t'i^,) (52) 

Jq dt'{Q\ exp f Zljggo ^i) ~ t')^i) ^t^i exp [t'Hq] \s) 

{ni)q = ^ ^ X -^^^ ^ • (53) 

Jo dt'{0\ exp (T^jego ^i) ^^P [^^i) 1^) 



B. Uncongested Regime 

We are now in a position to introduce the "uncongested" regime on a formal level. We will 
characterize this regime in terms of the existence of a special "universal" product-form solution 



to Eq. (48), as well as the finiteness of a particular expectation value (capturing the fact that no 
condensation of particles occurs at any nodes), at sufficiently large observational time t. Note that in 
the spirit of our operator-theoretic framework, we define the regime in terms of both the queueing 
network and the vector q on which one evaluates the generating function for occupation numbers 
in the network. Thus a given queueing network may be in the regime for some evaluations of its 
generating function (certain values of q) but not for others. This will then be related to belonging (or 
not belonging) to the regime for different types of large deviations (of current) through the standard 
Legandre transform, which maps the generating function (evaluated at different q) to the Cramer 
function (evaluated at different-sized deviations). 

We say that a given queueing network is in the "uncongested" regime for a given vector q if the 
ket vector \sq{n) is dominated by the ground state of Hq, i.e. at sufficiently large time, |sq(n)) ~ 
exp(— A(qr)t)|coh^(h.(qr)), holds. In other words, the spectrum of Hq is such that its ground state is 
separated from the excited states by a finite gap, and thus at the times much large than inverse value 
of the gap the solution is completely described by the ground state only, thus providing respective 
universality. We also require that (^i)q < oo for all nodes i of the network in the "uncongested" 
regime, ensuring that the expected number of particles (size of the queue averaged over time) does 
not diverge at any nodes as t — )■ oo. 
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The statistics of queueing networks in the "uncongested" regime may be analyzed using the tech- 
niques from Sectionjllj In particular, we substitute |Sq(n)) by ~ exp(— A(qr)t)|coh^(h.(qr)) in Eq. (48) 
to arrive at 



(54) 



yieQo: -hi{q) ^ A^^ + ^ qji\jihj{q) + \q, - \iQhi{q) = Q. (55) 



This set of relations generalizes the stationary (q = 1 ) relations ( 20pT ), and are thus consistent 
with A(l) = 0. Eqs. ( 54|[55 ) describe the right (ket) eigen-function of the ground-state of the 
evolution operator /Hamiltonian (49), while the corresponding left (bra)-eigenfunction is described 
in Appendix |A} 

Note that, replacing all internal g- variables by unity, the basic set of equations for h does not 
depend on the remaining goi and qio components. It follows that the respective h are identi- 
cal to the one derived before for the stationary ^(0) (no currents) setting. Moreover, A(qro) = 
Si {{qoi ~ l)Aoi — Koiqio — l)^j(O)). This observation translates (after the obvious Legandre trans- 
form) into the statement that all the currents entering and leaving the network are asymptotically 
Poisson, which (as already mentioned) was known previously in the Queuing literature [331 El] ■ 

For the sake of simplicity, hereafter we exclude the incoming and outgoing currents from consider- 
ation, which is achieved by setting qio = qoi = 1- The consistency of Eqs. (54) and Eqs. (55) (the two 
are just generalized versions of Eqs. ( 2lp0 )) translates into the following expression for the lowest 
eigenvalue: 



A(g) = - - 



(56) 



For sufficiently large t (it is our Large Deviation Parameter and thus should at least be significantly 
larger than the correlation time of the system) we obtain the following asymptotic expression for 



Pq{t) = ^(g)exp(-tA(qr)) ~ exp | t hi{q)\j{qij - 1) 

(i,i)eei 



*(q) = (0| exp 




|coh^(h(qr)). 



(57) 



(5^ 



Here, in evaluating the pre-exponential factor ^E'(g), we have used the assumption that the main 
contribution into Eq. (52) originates from t', t — t' ^ 1/A(qr). Thus, (0| exp(^j hi{q)di) in Eq. (|58|) 



is the bra- vector defined as the left eigenvector of Hg with the same eigenvalue A(qr) (see Appendix |A| 
for more details) ^. Here certain time- independent pre- factors (which do not impact the asymptotics 
up to exponential order) are ignored on the rhs. As we are interested in the statistics of the currents 



that scale (grow) linearly with t, one can replace the sum in Eq. (47) by an integral, invert the 



Note that expression for the analog of '!'(<?) eorrcspondent to Pq(n{t)) is significantly different; (0| exp ( X]iggo ^i) |cohm(fi.(q)) 
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relation, and arrive at the following saddle-point (large-deviation) expression: 



(i,i)egi 



\/{k,l)egiLk,lego: 3ki/qli= Yl 



dh,{q* 
dq, 



kl 



-{(fij - + hk{q*)\ku 



(59) 
(60) 



where S{j) is a convex function of its argument, also called the Cramer (or large-deviation) function. 



The dependence of the Cramer function on the current production j is defined implicitly via Eqs. (60 ) 



and Eqs. (55). 

Since A is fully defined by the solution of Eqs. (55), which is independent of m, the resulting 
expression for the Cramer function is also m-independent in the "uncongested" regime. This can- 
celation is quite remarkable. We note that the degeneracy is especially interesting, since it does 



not seem to extend to the time-independent pre- factor in P(J|t), ^{q) ( defined in Eq. (g ). 
A Queueing Theory interpretation is that in this regime the large deviations are not caused by 
the interactions of different particles in the network, and thus the same kind of deviations would 
have occured even if all nodes in the network were of M/M/oo type (as opposed to M/M/m with 
m < oo), in which different particles cannot interact and delay one- another in the network. ^. 

As explained above, the qualitative assumption that allowed us to extend the product-form ansantz 
to statistics of currents in Eqs. (54 60) was that the number of particles in the system (size of the 



queue averaged over time) did not diverge with time. For this assumption to hold, it must be the 
case that 



(0| exp (Ejeeo ^t^i rifeeeo 9mAhk{qW)\Q) 



< oo, 



0| exp ( EjGeo UkeGo 9m,ihkiq)a+)\0) 



(61) 



where h{q) and h{q) are solutions of Eqs. (p4, 55) and Eqs. (A2, A3) (which describe consistently 



the right /ket and left /bra ground state eigen-function of the evolution operator/Hamiltonian (49) 



We now comment on the phase transition which takes us from the uncongested to the congested 
regime. Picking a direction in the vector space of currents (dimensionality of the space is equal to 
the number of internal directed edges of the graph) and increasing the length of the vector in this 
direction, gradually (starting from the domain of values belonging to the "uncongested" regime) we 
will eventually reach the regime where congestion occurs, i.e. the condition (61) breaks down and 



particles accumulate over time at some nodes (condensation). To see this heuristically, we note that 
for large values of q, the main contribution to the generating function will be highly skewed towards 
very large values of current (due to the fact that q is raised to a power equal to the current), and the 
most likely way to attain such large values of currents will be for the queueing network to remain 
totally occupied over the entire time horizon, causing the number of particles in the network to 
diverge over time. Put in the context of large deviations, the most likely way to attain very large 



* Similar "mysterious" cancelation of the vorticity dependence was reported in '13' in the Cramer function of the entropy production for 
a polymer stretched by shear-vorticity flow. 
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values of currents is to have an accumulation over time of the number of particles in the network 
(condensation), leading to the system being overloaded over the entire time horizon, enabling all 
links in the network to generate current continuously over the entire time horizon. We focus in 
particular on the setting where the "breakdown" of the regime occurs initially at some particular 
unique node of the network, the expected number of particles in the system diverges at that node 
(over time), and a certain spectral condition holds at this breakdown point ( described in some detail 
below ). We will use these characteristics to formally define our "congested" regime. 

We now elaborate on the aforementioned spectral condition. In particular, the transition point 
from the "uncongested" to the "congested" regime has an interpretation in terms of the closing of 
the gap between the ground state and the low boundary of the continuous spectrum of Eq. (49). 
This "closing the gap" picture also assumes that other excited discrete states (which are present 
already in the case of a single station with feedback and m > 1, see the discussion in Section IV) 
do not cross the factorized ground state found above. This spectral interpretation also leads to 
an immediate conclusion/consequence in terms of the Cramer function shape at the values of the 
current that correspond to the considered uncongested-to-congested transition. Indeed, this "closing 
the gap" scenario translates into continuity of the Cramer function and its first derivatives at the 
transition. Stated differently (in the jargon of phase transition theory), the dynamical uncongested- 
to-congested transition with respect to the currents is second-order, and l/{ni)q plays the role of 
the order parameter at the congested node. 

Our final remark is about comparison of the distribution of queue averaged over time with re- 
spective distribution measured at the final moment t. The asymptotic analysis, described above for 
uncongested regime and extended in the following Subsection to congested regime, suggests that 
{ni{t))q is finite at any value of q and t — )■ oo, in particular for q = Qc where {ni)q = oo. Moreover, 
one conjectures that the two distribution functions, Pq{n) and Pq{n{t)), are different at any values 
of q except of the special one correspondent to the minimum of the Cramer function achieved at 
J = (J). This statement may be interpreted as a breakdown of ergodicity for any current but the 
special one correspondent to the steady distribution. 




FIG. 2: Transformation of the "overloaded" node, discussed in the text. The component of the graph, associated with the 
node, before and after the transformation are shown in the left-upper and low-right corners respectively. Dashed node on 
the modified graph correspond to new injections and departures with the Poisson rates exactly equivalent to the respective 
transition rates on the original graph. 
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C. Congested Regime 

We now consider a heuristic decomposition approach to extend the product-form description, 



utihzed above in the uncongested regime, beyond the domain bounded by Eqs. (61). We note that 
similar decomposition techniques have been considered throughout the Queueing theory hterature 
to understand how networks become congested |60l - [66l l68l 169] . 

Consider, for example, crossing the boundaries of Eq. (|6T| along a direction in the g-space, and 
thus violating the condition at some q^, at a single node, say i. Then, exactly at the point of crossing 
we can simply assume that the node is always congested (has an infinite queue in the waiting room), 
and thus the tellers at the node stay busy over the entire time horizon. This translates into the 
following obvious modification of the network graph: remove the node, and associate all the in/out 
edges for the node of the old graph to the new open-system in/out nodes with exactly the same rates. 
The transformation is shown in Fig. In the new domain the Eqs. (54 55) should be considered 
on the modified graph, and one needs to add the appropriate constant to A{q) from Eq. (56) to 
guarantee its continuity at q^,. Equivalently, one must compensate all probabilities to reflect the rare 
event that the given node remained busy over the entire duration. 

The general construction is illustrated in Fig. (|3]), and will also be illustrated in the next Subsection 
on our enabling example of a single node system with feedback. We note, of course, that in general 
the most likely way for rare events to occur in queueing networks may be quite complicated (see for 
example [22]), and what was explained above should be considered as an approximation which is 
not proven rigorously. 

We conclude this general construction by explaining the formal status of our derived results. All 
the asymptotic derivations so far, discussed in both the congested and uncongested cases, have 
relied on several important assumptions. The most important of these is the equivalence of the 
coherent state solution to the lowest eigen- value (ground state) solution. We have assumed that this 
(physically very plausible) assumption holds, and focused primarily on constructing the coherent 
state solution explicitly. We note, however, that our derivations have generally been non-rigorous 
in nature, and this assumption was not justified in any particular case. A considerable difficulty 
associated with formalizing our results (and justifying this assumption) lie in the fact that the 
coherent state approach does not allow us to analyze the entire spectrum of the operator. In view 
of the above, we find it important to validate this assumption, at least for a special case. We now 
proceed along these lines by studying the aforementioned case of a single station with feedback. 



D. Single station with feedback 



Perhaps the simplest example of a (non-trivial) queueing network (as discussed in Section IC) 



consists of a single station, one incoming channel/edge, one outgoing channel/edge, and one self- 
loop, all characterized by Poisson processes with rates A, fip and fi{l —p) respectively. The network 
is shown in the upper left corner of Fig. (|3|. 
In this example, the only interesting current is associated with the self-loop (as the others are 



Poissonian in the steady-state, see Section IC). We thus focus exclusively on the analysis of the 



current along the feedback arc, and the associated generating function (evaluated for the scalar q 
and corresponding current j). 
We start by considering the "uncongested" regime. Then Eqs. ( 54p6 ) become 



A 



/i(l - il-p)q)' 



A 



A(g-l)(l-p) 
l-{l-p)q 



(62) 
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FIG. 3: Original graph and transformed graph correspondent to uncongested and congested regimes for example of the single 
station with the feedback. 



According to Eqs. (59, 60), this results in 



3<3c- S{j) = X{l-p/2) ~ ^ +jln( 2j{l-p) 



To identify the "uncongested" regime breakdown we calculate h, according to Eq. (A3): 

p 



h 



1 - {l-p)q' 



Substituting into Eqs. (58) results in 



^(g) = 9m{Hq)h{q)), {n), 



d gm{h{q)h{q)z) 



dz gm{h{q)h{q)) 



(63) 



(64) 



(65) 



2 = 1 



From the above and Eq. (38), we find that the expected queue length remains bounded over time iff 



hh < m, and thus the critical value of congestion at which point one shifts regimes is 

1 

qc = 



1 — p 



(66) 




FIG. 4: Cramer function, S{j), of the feedback current shown for X = m — 1, p — 1/2 and fi — 3. Thick red and thick blue 
curves describe the uncongested and congested domains respectively. Dashed black line marks the value of jc- 
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The congested regime that occurs at q > Qc and j > jc corresponds to the situation when the single 
server is always occupied, and thus (conditioning on this event) jobs continually feedback over time 
as a m/x(l — p) standard Poisson process. Therefore, in the congested regime 



A 



— p)m{q — 1) + ( V A — ^Jp^m 



(67) 



Here the A — y/pjj/m 
and the condition of A-continuity at q 



term (a constant with respect to q) is found in accordance with Eq. (62) 



qc- Alternatively, in the congested regime A can 
identified with the lower edge s_ of the continuous spectrum of Hq, given by Eq. (74) . T he spectrum 
of Hq for a simple model under consideration is analyzed in some detail in Section IV Performing 



the Legendre transform on Eq. (67), we arrive at the following expression for the Cramer function 



of the feedback current in the congested regime 



J > jc 



S{j) = (^y/X — ^/pjJm^ + m/i(l — p) + j ln( 



em/i(l — p) 



{61 



Comparing Eq. (63) with Eq. (68) we also observe that the Cramer function is smooth (first derivative 
is continuous) at j = jc, thus confirming that the dynamical phase transition described here is of the 
second order (continuous). The change of the Cramer function shape across the transition is shown 
in Fig. m 

The obtained expression for the Cramer function in the congested regime Eq. (68) has a very 
simple and transparent interpretation. Since the server is always occupied, the feedback current 
generation is a standard Poisson process, so that S{j) = fim{l — p) + j ln(j'/(em/i(l —p))), where 
g-fcSo jg ^]-^g probability to keep the server busy, presented with exponential accuracy. The latter 



is given by the probability e 



of creating the incoming and outgoing currents of the same 



value (x! to provide the marginal stability of the server, maximized with respect to to. Since both 
currents are generated by independent Poisson processes with the rates A and pfim, respectively, 
we have So{u)) = A + pmfi + u ln(aj^/ {Xpfim)e'^), and minimization with respect to u results in 



So = {y/X — ^ppmY, which reproduces Eq. (68) 



IV. DIRECT ANALYSIS OF THE SINGE-STATION FEEDBACK SYSTEM 



As shown in the previous Sections, an understanding of the evolution operator ground state struc- 
ture is essential for predicting long-time statistics of currents. We have also argued that the ground 
state can be described analytically for networks with general graphical structure. However, these 
arguments were indirectly dependent on certain information about the rest of the spectrum - specif- 
ically on the fact that the ground state is separated from the continuous spectrum by a gap which 
collapses at the dynamical phase transition. In general, information on the entire spectrum is dif- 
ficult to obtain. The main point of this Section is to gain a broader understanding of the simple 
single-node network with feedback along these lines. Therefore in this Subsection we analyze the 
full spectrum of the problem with the single feedback and we confirm the general picture of the 'gap 
emergence', and collapse at the phase transition point suggested above on the basis of only partial 
(ground state) analysis. 

Our starting point is the dynamical equation Pq{n;t) = Yl%.o1^ P^^i j'^'^) ^^r the generating func- 
tion of current j via the feedback arc. Here P{n,j;t) is the joint probability distribution function 
of the number of particles, n standing at the queue at the time t, and the number of particles, j, 
that have passed through the feedback loop by time t. Following directly the proper generalization. 
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FIG. 5: Spectrum of the evolution operator as a function of q for m = 5, A = 0.5, /i = 2.0, p = 0.5 and A'' = 50. The lower 
bottom line corresponds to the ground state solution A(g). 



according to Eqs. (43 49) of the ME (33) we derive: 
d 



dt 



P,in; t) = X {Pg{n - 1; t) - Pg(n; t)) + ^ip {Omin + l)P,(n + 1; t) - e„(n)Pg(n; t)) 



+fi{l - p){q - l)eUn)Pg{n;t), 



(69) 



where the last term on the r.h.s. accounts for the current. Performing the Laplace transform of 
Pq{n;t) over time, Ps-q{n) = exp{st)Pq{n;t)dt (with s considered as a spectral parameter), we 
arrive at the following spectral equation: 



A - fxpOmin) + -p){q- l)6min))P,.g{n) 



-\Ps-q{n 



I)- ^pe^{n + l)P,,q{n + l). (70) 



Here the inhomogeneous part, dependent on the initial condition at t = 0, is ignored. 
The "boundary" conditions over n for Eq. (70) read P< 



s\q\ 



and P.-nin — > oo) 



0. We 



will relax the last condition, assuming the following natural finite waiting room regularization: 
Ps:q{N + m) = for some sufficiently large value of N. We will study the spectrum at finite A^, and 
show towards the end of the calculations that the iV — )■ oo limit is well defined. We note that the 
particular choice of regularization turns out to be unessential for the limits we consider. In order to 



find the eigenvalues (spectrum) we solve the recurrence relation (70) for general values of s, and use 
the aforementioned boundary condition, and the normalization condition Ps-q{0) = 1. To summarize, 
the set of conditions that complement Eq. ((70| to provide the full description of the spectrum are 



0, Ps;,(0) = l, Ps-q{N + m) 



0, 



(71) 



and we are mainly interested in studying the N ^ oo limit. 

A numerical solution of the eigenvalue problem is illustrated in Fig. ([s]) for m = 5, A = 0.5, fi = 
2.0, p = 0.5 and = 50. The simulations show emergence at g = 1 of three discrete eigenstates well 
separated from the (still discrete) band of states. We tested numerically as one increases A^, the 
q = value of the ground state s (top line in Fig. |5]) approaches 0. We observed that the asymptotic 
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dependence of the ground state eigen- value on s is fully consistent with the prediction of the ground- 



state coherent-state theory explained above in Section HID, specifically with Eq. (62). Moreover 



our prediction of Qc (where the gap between the ground state and the continuous band collapses), 



as given by Eq. (66), is fully consistent with the respective dependencies of the gap collapse in our 
computations. Let us also mention that the two non-ground state eigen-values observed at g = 1 
never cross with each other or with the ground state at g > 1, and merge into the continuous band 
(consequently one after another) at values of q smaller than Qc- Experimenting with larger m, we 
observe that the discrete spectrum becomes equidistant in the m — )■ oo limit, which is fully consistent 
with the interaction-free nature of the limit. 

We now analyze the spectrum analytically. Let us begin by making some preliminary observations. 
First, the linear dependence on s implies that the function Ps-q{n) is an n-th order polynomial in 



s, suggesting that there are always exactly N + m solutions of the last condition in Eq. (71) 



Ps:q{N + m) = 0. Second, at n > m the factor 6m{n) becomes constant and the recursive relations 



(70) can be solved analytically. Thus, looking for solution of Eqs. (70) at n > m in the Ps-q{n 



c+p^ + c^p"i form, one obtains the following expression for p± (which depends on g, s and the rates): 

p± = (\ + ppm — p(l — p)(q — l)m — s ± a/(A + ppm — p(l — p)(q — l)m — sV — iXpprn] . 

zppm \ / 

(72) 

The boundary condition Ps-q{N + m) = translates into the following condition ( dependent on 
Ps-g{m) and Ps;q{rn - 1)): 

I ^) _ (Ps-Am) - p.Ps.qjm - l))p^ + {P,,q{m) - p^Ps.qjm - l))pj! _ ^ 

'"^ p+- p- 

The above has different types of solutions dependent on the values of s with respect to the following 
two threshold values: 



s± = A + ppm + yu(l — p){q — l)ra ± ^/AXppm = (^^/X ± y/p^pm^ — p(l — p){q — l)m. (74) 

At s < S-, both eigenvalues p± are real and p+ > p_. In this case, as we are interested in the 
N ^ oo limit, one can simply ignore the p^ contributions in (73), thus leading to the following 
relation: Ps-^q{m) — p^Ps-qlm — 1) = 0. This then replaces the last condition in Eq. (71), w.r.t. 
describing Ps;q{n) at < n < m. We were not able to solve this reduced system of equations 
analytically at any values of m, and thus to find the spectrum in its full glory. However, the 
information just provided is already sufficient for a heuristic derivation of the lowest eigenvalue, and 
the value of the gap between it and the continuous band. 
We focus on a particular single form solution of the reduced system of equations for Ps-^q{n) with 



< n < m, corresponding to s = A and Ps-q{n) = /n\ for n < m, with h and A from Eq. (62). In 
this case, one has p+ = A/ (pph) and p_ = h/m, so this is indeed the ground state solution. However, 
let us note for the sake of accurateness that we have not formally proven that this special solution 
always corresponds to the lowest eigenvalue (the ground state). However, this is exactly what we 
observed in our numerical experiments with different values of m. 
Eq. ((73j) also allows the continuous spectrum of the operator to be identified. At s_ < s < s+ we 



have |p+| = |p_|, and one has to keep both terms in Eq. (73). As long as \Ps-q{m + 1) — p^Ps-q{m 



|-Ps;g('^) ~ P+Ps;q{.'m — 1)|, Eq. (73) has at least solutions. At A^ ^ 1 this region turns into the 
continuous band of the spectrum. We conjecture that the system of equations does not have any 
solutions for s > s_|_. (Once again, this is confirmed in simulations but we do not have an explicit 
way of proving it.) 
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Finally, we conclude that the transition between the uncongested regime and congested regime 



takes place at A = s_ which corresponds exactly to q = qc ( from Eq. (66) ). As one can easily see 



the value of h at this point is equal to h = \f\fijup < 1, so the queue length in the final moment 
does not diverge. 

V. CONCLUSIONS AND PATH FORWARD 

Let us briefly recall the highlights of this manuscript. Our analysis was focused on the the gen- 
erating function of currents over a Jackson (queueing) network, with an eye towards analyzing how 
large currents accumulate over time. We began by giving some relevant background in Queueing 
Theory, and discussing some tie-ins with recent work in statistical physics. We then adopted the 
Doi-Peliti technique for describing the dynamics of the Jackson network. We used this formalism to 
show that the ground state of the respective evolution operator has a well-defined and analytically 
tractable pro duct /coherent state form in a particular regime. These results were translated into 
an implicit analytical expression for the Cramer function of currents in this "uncongested" regime, 
where the ground state was well separated by a gap from the continuous spectrum. We also observed 
that crossing the surface in the phase space of observed currents, where the spectral gap collapses, 
corresponds to the congestion of a node in the network. We suggested a heuristic graph-reduction 
scheme which allows the statistics of currents to be described in this partially congested setting. 
Finally, we validated the general results using the example of a single node feedback system, where 
many of our assumptions could be verified directly by performing an explicit analysis of the evolution 
operator spectrum. 

We consider this study more like an opening for further exciting research along the following lines: 



Following the discussion in Section [III C[ we naturally conjecture that in a large network, 
gradually increasing the observed current (s) will lead to a transition from the uncongested 
to the fully congested regime via a number of steps, each characterized by an increase in 
the number of congested nodes. It would be interesting and important to explore the specific 
sequence of phase transitions separating the space of observed currents into cells. A particularly 
interesting question concerns the algorithmic complexity of identifying these cells and exploring 
their geometric structure (e.g., possible convexity). 

We have considered a queuing model with the transition rates independent of the node occupa- 
tion numbers. This condition can be relaxed and such an extension of our theory should allow 
for non-uniform dependence of the rates on the occupation numbers. 

Extensions of the current-statistics theory to the so-called multi-class networks, where different 
classes of particles are treated differently at the stations (for example having different priorities) 
are less trivial, yet we anticipate such a generalization is still possible. 

In this manuscript we have considered solely open queuing networks. It would be interesting 
and instructive to generalize the theory of current statistics to the case of closed (particles are 
not leaving or entering the network) and semi-open (particles are injected into the system and 
leaving it, but in such a way that the total number of particles is conserved) networks. 

We may also consider networks of fixed structure, yet with rates changing in time, for example 
in a periodic fashion. To describe statistics of currents in this case, and especially in the 
regime where the typical correlation time of the rate changes are comparable to the inverse 
rates, constitutes another interesting future challenge. (Note that an approach blending the 
techniques described in this manuscript with the ones discussed recently in the context of the 
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so-called Jarzynski equality [70] and work relations [7T1 [72], both closely related to the subject 
of fluctuation theorems p^l - [T2] . may prove fruitful for this task.) 

• The assumption of infinite waiting room was crucially important for advancing the prod- 
uct/coherent state decomposition to the statistics of currents. Analyzing the current statistics 
in the regimes when all or some waiting rooms are of a finite capacity is yet another challenging 
task, as this finiteness brings in a new type of inter-particle interaction. 

• It would be an interesting challenge to put the ideas presented in this paper on a more rigorous 
mathematical foundation. This would aid greatly in understanding the formal relationship 
between the regimes identified in this paper and the sample path large deviations properties 
of queueing networks. We note that the authors are currently undertaking preliminary work 
along these lines, with the work directed more towards the Queueing Theory community. 

• Our analysis suggests that in large queuing networks one may expect emergence of multiple 
transitions as one increases the current. More generally, and in the spirit of [TSj [71], it would 
be interesting to explore the field of the so-called qualitative queuing network theory via the 
methods/techniques discussed in this manuscript. 

• Finally, all of the above should be used not only to study existing (man- or nature- made) 
networks, but also to guide construction of future technological networks with desired properties. 
In other words, we suggest to use this analysis for control and optimization of networks in new 
areas such as power, and even more generally energy, distribution. 



VI. ACKNOWLEDGMENTS 



We are thankful to David Gamarnik for consulting us on many issues related to Queuing Theory, 
and Sergey Foss, Bill Massey and Alexander Rybko for enlightening conversations. This mate- 
rial is based upon work supported by the National Science Foundation under CHE-0808910 (VC) 
and CCF-0829945 (MC via NMC). The work at LANL was carried out under the auspices of the 
National Nuclear Security Administration of the U.S. DoE at LANL under Contract No. DE-AC52- 
06NA25396. KT acknowledges support of an Oppenheimer Fellowship at LANL, and DAG work on 
the project was a part of his summer internship (GRA program) at LANL. 



Appendix A: Left Ground State of the Hamiltonian 



In this Appendix we construct the left ground state of the Hamiltonian (49). We start by recalling 
that according to Eq. ([7|, the bra-vector (0| exp(^^ Sj) is the left zero eigen-function of the Hamil- 
tonian (49) at q = 1. However, it ceases to be an eigen- vector at g 7^ 1. On the other hand, it 



is easy to check that the "exponential" bra- vector (0|exp(/ia) is in fact a left eigen- vector of the 
creation operator a+, with the eigen-value h, (0| ex-p{ha)a~^ = h{0\ exp{ha). This suggests searching 
for a left eigen- vector of the Hamiltonian (49), in the exponential form: 



(0|exp ^hitti 



(Al) 
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Then, from {sq\Hq = —A{q){sq\ and Eq. (49) (combined with utihzing the aforementioned feature 
of the creation operators) one derives the following set of conditions on the h vector (of c- numbers): 



"^{qoi - hi{q))Xoi = A(g) 



(i,i)6Gi 



: Xio {qio ~ hi{q)) + ^ijiQijhjiq) - hi{q)) = 0. 



(A2) 
(A3) 



It is straightforward to verify that the resulting A{q) and A(qr) from Eq. (54) are identical, A(qr) 
A(g). 
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